Análisis de Series Temporales de los Espectadores en las salas de cine
Alejandro Vaca Serrano
12 de marzo de 2019
Introducción: Análisis previo de los datos
En este apartado nos centraremos en entender la composición de nuestros datos, así como la tendencia, estacionariedad y estacionalidad de los mismos, con el fin de ajustar posteriormente diversos modelos predictivos para finalmente quedarnos con el más efectivo a la hora de predecir. Por el camino realizaremos algunos test estadísticos y gráficos que nos ayudarán a determinar cuál debe ser el tratamiento de la serie, prevaleciendo siempre el criterio de minimizar el error de predicción.
## Warning: package 'ggplot2' was built under R version 3.5.2
Como vemos en el gráfico, la tendencia es claramente descendente, mostrando que en los últimos años el uso del cine como entretenimiento ha decaído considerablemente. Una limitación de esta serie temporal es que llega sólo hasta el año 2012 (este año está completo), justo el año en el que entraron empresas como Netflix al mercado. Esto es una limitación en el sentido de que no podremos incluir esta variable (si existe o no netflix o hbo y el número de usuarios). Sin embargo sí que podremos acceder, no sin dificultad, a los datos del precio medio de una entrada del cine en España, en este link. En el gráfico superior podemos ver algunas estacionalidades, ya que los picos de mayor uso de las salas de cine son más o menos cíclicos. Lo comprobaremos a lo largo de la práctica observando las funciones de autocorrelación total y parcial, entre otros métodos que iremos viendo en el análisis.
Separaremos la muestra en datos de train y test, dejando como entrenamiento todos los datos previos al año 2009 y los 3 años restantes como período de prueba del modelo.
## Warning: package 'tseries' was built under R version 3.5.2
En los gráficos superiores podemos observar la descomposición de la serie en los diferentes elementos que la conforman. La tendencia en el tiempo es claramente descendente, mostrando la progresivamente menor popularidad de este medio de ocio. Podemos observar, además, que existe una cierta estacionalidad en los datos, pues los mismos “picos” de usuarios se repiten a lo largo del tiempo; el más grande de ellos se repite cada año. Es doble, pues hay al rededor de dos picos por año, con una bajada también grande antes de estos dos picos sucesivos. Respecto al ruido, parece oscilar al rededor de una media, aunque no es del todo constante en el tiempo en varianza, ya que hacia el año 2008 la varianza comienza a aumentar. Esto es positivo a la hora de aplicar un modelo GARCH, pues cuando las series tienen una homocedasticidad marcada este tipo de modelos no tienen tanto sentido, ya que la varianza del mes siguiente va a ser prácticamente igual a la del mes actual.
library(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recodelibrary(MASS)## Warning: package 'MASS' was built under R version 3.5.2##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## selectbox_cox <- boxcox(users ~ date,
data = train,
lambda = c(0, 0.5, 1))
lambda <- box_cox$x[which.max(box_cox$y)]
lambda## [1] 0.6666667powerTransform(train.ts)## Estimated transformation parameter
## train.ts
## 0.6566137Vemos que \(\lambda\) está entre 0.5 y 1, por lo que tanto diferenciar como hacer la raíz cuadrada podrían ser buenas ideas. Lo que haremos más adelante será probar todo y quedarnos con el modelo que nos dé un menor error en la predicción, centrándonos más en el performance predictivo que en la robustness estadística del modelo, pues para ponerlo en producción lo que más dinero ahorrará o generará para la empresa será predecir de la mejor manera posible.
Funciones de Autocorrelación Parcial y Total: Ajuste de p y q.
Recordemos que la serie es mensual, por lo tanto el 1 del retardo hace referencia a un retardo de 12 meses (el eje X debemos multiplicarlo x 12 si queremos tenerlo en meses). Vemos en el gráfico de correlación parcial de la serie que existe una correlación significativa con el mes anterior, con 6 meses atrás, con 11 meses atrás y con 12 meses atrás. Esto haría referencia a la q del modelo, es decir la parte de la media móvil. Por otro lado vemos en la ACF que existe correlación con los datos de 1 año, 2 años y 3 años antes. También parece haber correlación con 1 mes antes, 8, 11 y 12 meses antes. Como vemos este comportamiento de que se correlacionan las observaciones con su homólogo de 1 año, 2 y 3 años antes, vamos a incluir una parte estacional.
Probamos primero a ajustar la parte estacional del modelo, de orden 1.
Primer Ajuste
## Series: train.ts
## ARIMA(0,0,0)(1,0,0)[12] with non-zero mean
##
## Coefficients:
## sar1 mean
## 0.7340 10760.4226
## s.e. 0.0713 452.7668
##
## sigma^2 estimated as 2397087: log likelihood=-844.96
## AIC=1695.92 AICc=1696.18 BIC=1703.61## sar1 intercept
## sar1 1.00000000 -0.04958451
## intercept -0.04958451 1.00000000## Retard p-value
## [1,] 6 0.10506178
## [2,] 12 0.05530046
## [3,] 18 0.07325149
## [4,] 24 0.13259838
## [5,] 30 0.06617838
## [6,] 36 0.04317806
## [7,] 42 0.10556828
## [8,] 48 0.08014467
Vemos que se nos salen aún algunas correlaciones parciales de las bandas (de Barlett); nos llama especialmente en el PACF la banda del 1, ya que esta se sale mucho. En el caso de la ACF también se sale bastante esta misma banda. Como vemos en la PACF que existe correlación significativa con el valor de la serie en el mes anterior, decidimos probar un ARMA(0, 0, 1)x(0, 0, 1). Algo positivo es que el p-value de las correlaciones nos sale significativo para casi todos los lags, con lo cual estamos más cerca de cumplir con los supuestos del modelo ARMA.
Segundo Ajuste
ajuste2 <- Arima(train.ts,
order = c(0,0,1),
seasonal = list(order = c(1,0,0), period = 12),
method = "ML")
ajuste2## Series: train.ts
## ARIMA(0,0,1)(1,0,0)[12] with non-zero mean
##
## Coefficients:
## ma1 sar1 mean
## 0.3080 0.6986 10761.9363
## s.e. 0.1008 0.0770 518.2067
##
## sigma^2 estimated as 2233367: log likelihood=-840.47
## AIC=1688.95 AICc=1689.39 BIC=1699.21## ma1 sar1 intercept
## ma1 1.000000000 -0.20087860 0.002181477
## sar1 -0.200878601 1.00000000 -0.047051282
## intercept 0.002181477 -0.04705128 1.000000000## Retard p-value
## [1,] 6 0.7389242
## [2,] 12 0.2912638
## [3,] 18 0.5169280
## [4,] 24 0.7311593
## [5,] 30 0.4436429
## [6,] 36 0.3183765
## [7,] 42 0.4317232
## [8,] 48 0.4415686
Vemos que, especialmente en la PACF, persiste la correlación significativa con el valor del año anterior, por lo tanto decidimos incluir un nuevo componente estacional, y probaremos un ARMA(0, 0, 1)x(1, 0, 1).
Tercer Ajuste
ajuste3 <- Arima(train.ts,
order = c(0,0,1),
seasonal = list(order = c(1,0,1), period = 12),
method = "ML")
ajuste3## Series: train.ts
## ARIMA(0,0,1)(1,0,1)[12] with non-zero mean
##
## Coefficients:
## ma1 sar1 sma1 mean
## 0.3366 0.9012 -0.4260 10748.0593
## s.e. 0.1029 0.0492 0.1183 621.2377
##
## sigma^2 estimated as 1984567: log likelihood=-835.7
## AIC=1681.4 AICc=1682.07 BIC=1694.22## ma1 sar1 sma1 intercept
## ma1 1.00000000 0.08343816 -0.25552459 0.00475227
## sar1 0.08343816 1.00000000 -0.68785287 -0.00188364
## sma1 -0.25552459 -0.68785287 1.00000000 -0.02603497
## intercept 0.00475227 -0.00188364 -0.02603497 1.00000000## Retard p-value
## [1,] 6 0.2014748
## [2,] 12 0.1403626
## [3,] 18 0.2180477
## [4,] 24 0.4550048
## [5,] 30 0.3411850
## [6,] 36 0.2758098
## [7,] 42 0.4415517
## [8,] 48 0.4268327
Vemos que en la PACF se sigue saliendo una barra de los límites, en concreto la correlación de la observación en el período t con el valor de 2 años y 6 meses atrás… Podríamos probar algunos cambios en los órdenes del modelo para ver si esto desaparece, aunque no hay ninguno a primera vista que pudiera ser de ayuda. Una idea sería meterle una diferencia estacional. Una buena noticia es que seguimos sin tener que diferencir la serie obligatoriamente, pues los intervalos de confianza de nuestros coeficientes no incluyen al 1.
Cuarto Ajuste
Después de realizar muchas pruebas, nos encontramos con que, después del ARMA anterior que mostramos, el siguiente mejor sería un arima: ARMA(0, 1, 1)x(0, 1, 1). Sin embargo este modelo no mejoraría esto que comentábamos de que sigue existiendo una correlación parcial con las observaciones de dos años y medio atrás. Más adelante meteremos efectos de calendario y similares, y esperamos que estos nos puedan ayudar a explicar este suceso.
ajuste4 <- Arima(train.ts,
order = c(0,1,1),
seasonal = list(order = c(0,1,1), period = 12),
method = "ML")
ajuste4## Series: train.ts
## ARIMA(0,1,1)(0,1,1)[12]
##
## Coefficients:
## ma1 sma1
## -0.9185 -0.7404
## s.e. 0.0630 0.1364
##
## sigma^2 estimated as 1563745: log likelihood=-714.6
## AIC=1435.2 AICc=1435.51 BIC=1442.46## ma1 sma1
## ma1 1.00000000 -0.09393753
## sma1 -0.09393753 1.00000000## Retard p-value
## [1,] 6 0.5054048
## [2,] 12 0.6491118
## [3,] 18 0.4477164
## [4,] 24 0.6099739
## [5,] 30 0.6912063
## [6,] 36 0.7900751
## [7,] 42 0.8462380
## [8,] 48 0.9300662
Probando con la raíz cuadrada:
Primer Ajuste
ajuste1_sqrt <- Arima(train2.ts,
order = c(0,0,0),
seasonal = list(order = c(1,0,0), period = 12),
method = "ML")
ajuste1_sqrt## Series: train2.ts
## ARIMA(0,0,0)(1,0,0)[12] with non-zero mean
##
## Coefficients:
## sar1 mean
## 0.7353 103.0976
## s.e. 0.0718 2.2057
##
## sigma^2 estimated as 56.43: log likelihood=-333.46
## AIC=672.92 AICc=673.18 BIC=680.61## sar1 intercept
## sar1 1.000000 -0.061306
## intercept -0.061306 1.000000## Retard p-value
## [1,] 6 0.11317365
## [2,] 12 0.07335951
## [3,] 18 0.09586430
## [4,] 24 0.16246055
## [5,] 30 0.09605351
## [6,] 36 0.06446409
## [7,] 42 0.15921812
## [8,] 48 0.11712114
Segundo Ajuste
## Series: train2.ts
## ARIMA(1,0,0)(1,0,0)[12] with non-zero mean
##
## Coefficients:
## ar1 sar1 mean
## 0.3102 0.7042 103.1031
## s.e. 0.0982 0.0758 2.8135
##
## sigma^2 estimated as 52.25: log likelihood=-328.74
## AIC=665.48 AICc=665.92 BIC=675.74## ar1 sar1 intercept
## ar1 1.000000000 -0.13106316 -0.001174168
## sar1 -0.131063164 1.00000000 -0.051907019
## intercept -0.001174168 -0.05190702 1.000000000## Retard p-value
## [1,] 6 0.8491386
## [2,] 12 0.3958939
## [3,] 18 0.6000988
## [4,] 24 0.7914190
## [5,] 30 0.4284576
## [6,] 36 0.3247987
## [7,] 42 0.4542195
## [8,] 48 0.4616925
Tercer Ajuste
## Series: train2.ts
## ARIMA(1,0,0)(1,1,1)[12]
##
## Coefficients:
## ar1 sar1 sma1
## 0.3471 -0.3230 -0.2515
## s.e. 0.1157 0.1798 0.1898
##
## sigma^2 estimated as 45.24: log likelihood=-279.82
## AIC=567.65 AICc=568.15 BIC=577.37## ar1 sar1 sma1
## ar1 1.0000000 0.1754495 -0.4071704
## sar1 0.1754495 1.0000000 -0.7505946
## sma1 -0.4071704 -0.7505946 1.0000000## Retard p-value
## [1,] 6 0.11419890
## [2,] 12 0.13121151
## [3,] 18 0.09610893
## [4,] 24 0.25061455
## [5,] 30 0.29791174
## [6,] 36 0.45922481
## [7,] 42 0.63002688
## [8,] 48 0.59753012
Parece haber una pequeña anomalía, pues encontramos una correlación parcial con las observaciones de hace año y medio significativa; sin embargo todo el resto de “barras” (correlaciones) se encuentran dentro de las bandas de Barlett.
Elevando la Serie a Lambda (Approximate Box-Cox Transformation)
Primer Ajuste
ajuste1_lambda <- Arima(train3.ts,
order = c(0,0,0),
seasonal = list(order = c(0,1,0), period = 12),
method = "ML")
ajuste1_lambda## Series: train3.ts
## ARIMA(0,0,0)(0,1,0)[12]
##
## sigma^2 estimated as 1.111: log likelihood=-123.62
## AIC=249.23 AICc=249.28 BIC=251.66
Segundo Ajuste
ajuste2_lambda <- Arima(train3.ts,
order = c(0,0,0),
seasonal = list(order = c(0,1,1), period = 12),
method = "ML")
ajuste2_lambda## Series: train3.ts
## ARIMA(0,0,0)(0,1,1)[12]
##
## Coefficients:
## sma1
## -0.3652
## s.e. 0.0941
##
## sigma^2 estimated as 0.9496: log likelihood=-117.37
## AIC=238.74 AICc=238.89 BIC=243.6
Tercer Ajuste
Después de probar con varios, el modelo que mejor se ajusta es el SARIMA(0, 1, 1)x(1, 1, 1)[12]. Los residuos de este modelo presentan aún algunas anomalías, como veremos en los gráficos ACF y PACF, y es que no conseguimos que todas las correlaciones sean no significativas (dentro de la línea de Barlett) con \(y^\lambda, \lambda = 0.666667\). Esto podría ser debido a outliers, lo cual examinaremos más adelante.
ajuste3_lambda <- Arima(train3.ts,
order = c(0,1,1),
seasonal = list(order = c(1,1,1), period = 12),
method = "ML")
ajuste3_lambda## Series: train3.ts
## ARIMA(0,1,1)(1,1,1)[12]
##
## Coefficients:
## ma1 sar1 sma1
## -0.9043 -0.3177 -0.5540
## s.e. 0.0554 0.1595 0.1745
##
## sigma^2 estimated as 0.6321: log likelihood=-103.11
## AIC=214.21 AICc=214.72 BIC=223.89
Análisis de Intervenciones: Variables Explicativas.
Ahora vamos a añadirle a los mejores modelos de cada tipo (de cada transformación que hemos probado) las variables exógenas. Éstas van a ser relativas a si cae en semana santa, al número de días festivos del mes, al precio del cine medio en España ese año, y la variable dt, que se puede ver definida en el código.
En el chunk de código inferior vamos a definir una función que nos de las variables explicativas. Aparte de los festivos y demás efectos de calendario comunes a casi todas las series temporales, le añadiremos algunas concretas de este caso. Hemos cogido como “marco” la función definida por Daniel Vélez, y le hemos realizado diversos cambios para incorporar, por ejemplo, el efecto del día 2 de Mayo. Nos gustaría haber incluido que el día fuera miércoles como efecto, pues estos días hay descuento en los cines, pero es desde el año 2014 así que nos olvidamos. Sí vamos a incluir el precio medio del cine en cada año, pues esto esperamos que tenga una relación directa con el número de espectadores en las salas de España. Sacar estos datos nos llevará algo de trabajo, pues el único estudio disponible es este y no tiene una tabla de la que podamos sacar los datos, sino que estos tienen que ser sacados a mano. De igual manera, los introducimos como variables explicativas.
Primer Modelo
xregs = as.matrix(expl_train[ , c("semanaSanta", "dt", "bisiesto", "festivo", "price")])
ajuste5 <- Arima(train.ts,
order = c(0,0,1),
seasonal = list(order = c(1,0,1), period = 12),
xreg = xregs,
method = "ML")
ajuste5## Series: train.ts
## Regression with ARIMA(0,0,1)(1,0,1)[12] errors
##
## Coefficients:
## ma1 sar1 sma1 intercept semanaSanta dt bisiesto
## 0.4353 0.9492 -0.6054 13895.825 360.7003 -113.3265 333.2582
## s.e. 0.0846 0.0362 0.1313 1313.897 74.9046 29.8079 690.8709
## festivo price
## 909.0092 -737.4628
## s.e. 350.3033 208.0835
##
## sigma^2 estimated as 1270283: log likelihood=-812.58
## AIC=1645.16 AICc=1647.75 BIC=1670.8summary(ajuste5)## Series: train.ts
## Regression with ARIMA(0,0,1)(1,0,1)[12] errors
##
## Coefficients:
## ma1 sar1 sma1 intercept semanaSanta dt bisiesto
## 0.4353 0.9492 -0.6054 13895.825 360.7003 -113.3265 333.2582
## s.e. 0.0846 0.0362 0.1313 1313.897 74.9046 29.8079 690.8709
## festivo price
## 909.0092 -737.4628
## s.e. 350.3033 208.0835
##
## sigma^2 estimated as 1270283: log likelihood=-812.58
## AIC=1645.16 AICc=1647.75 BIC=1670.8
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -97.86562 1072.937 849.0645 -2.007611 8.058349 0.650576
## ACF1
## Training set 0.06290485coeftest(ajuste5)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 4.3529e-01 8.4565e-02 5.1474 2.641e-07 ***
## sar1 9.4918e-01 3.6179e-02 26.2354 < 2.2e-16 ***
## sma1 -6.0538e-01 1.3127e-01 -4.6118 3.992e-06 ***
## intercept 1.3896e+04 1.3139e+03 10.5760 < 2.2e-16 ***
## semanaSanta 3.6070e+02 7.4905e+01 4.8155 1.469e-06 ***
## dt -1.1333e+02 2.9808e+01 -3.8019 0.0001436 ***
## bisiesto 3.3326e+02 6.9087e+02 0.4824 0.6295403
## festivo 9.0901e+02 3.5030e+02 2.5949 0.0094613 **
## price -7.3746e+02 2.0808e+02 -3.5441 0.0003940 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Box.test.2(residuals(ajuste5),
nlag = c(6,12,18,24,30,36,42,48),
type="Ljung-Box")## Retard p-value
## [1,] 6 0.7827383
## [2,] 12 0.8363230
## [3,] 18 0.7038717
## [4,] 24 0.8862800
## [5,] 30 0.9049797
## [6,] 36 0.8893456
## [7,] 42 0.8556661
## [8,] 48 0.9064435cor.arma(ajuste5)## ma1 sar1 sma1 intercept semanaSanta
## ma1 1.00000000 0.03887024 -0.094470090 -0.04136217 -0.029139960
## sar1 0.03887024 1.00000000 -0.825996340 0.31795618 0.014912022
## sma1 -0.09447009 -0.82599634 1.000000000 -0.41203364 0.001940848
## intercept -0.04136217 0.31795618 -0.412033638 1.00000000 -0.031848228
## semanaSanta -0.02913996 0.01491202 0.001940848 -0.03184823 1.000000000
## dt 0.07465806 -0.02273704 0.038908550 -0.01923734 0.142908868
## bisiesto 0.07091865 0.04288069 -0.094600353 0.06675116 0.055498652
## festivo -0.01469069 0.04987521 -0.068116890 -0.19803690 0.058844929
## price 0.05357171 -0.34904465 0.450341177 -0.88010265 -0.008396677
## dt bisiesto festivo price
## ma1 0.07465806 0.07091865 -0.01469069 0.053571711
## sar1 -0.02273704 0.04288069 0.04987521 -0.349044653
## sma1 0.03890855 -0.09460035 -0.06811689 0.450341177
## intercept -0.01923734 0.06675116 -0.19803690 -0.880102650
## semanaSanta 0.14290887 0.05549865 0.05884493 -0.008396677
## dt 1.00000000 -0.01162851 -0.03982757 0.024539631
## bisiesto -0.01162851 1.00000000 0.07225523 -0.107366327
## festivo -0.03982757 0.07225523 1.00000000 -0.053487506
## price 0.02453963 -0.10736633 -0.05348751 1.000000000Tenemos alguna correlación algo alta, como entre el precio y el intercept (-0.88); de todas formas la dejaremos de momento. Respecto a las demás variables, parece que la única que no es estadísticamente significativa es la binaria bisiesto. Una mala noticia es que, como vemos, el AIC baja poco con respecto al mismo modelo ARMA pero sin regresores externos. Vamos a volver a dibujar el ACF y PACF de los residuos del modelo, para comprobar si hemos conseguido resolver esa pequeña anomalía que observábamos (la correlación significativa con t - 2.5 años) (t - 30 meses, pues la serie es mensual). Antes de hacer esto debemos eliminar la variable no significativa (bisiesto).
Parece que sí ha quedado resuelto el pequeño conflicto que existía. Ahora ya tenemos todas las correlaciones dentro de las barras de Barlett y por lo tanto, teóricamente hemos capturado la información existente en la serie temporal.
Segundo Modelo: Con la Raíz Cuadrada
Como nos da problemas al sacar la diferencia estacional, cambiamos ligeramente el modelo presentado anteriormente con la Raíz Cuadrada. Usaremos un SARIMAX(1,0,0)x(1,0,1).
ajusteSqrt = Arima(train2.ts,
order=c(1,0,0),
seasonal=list(order=c(1,0,1), period=12),
xreg=xregs,
method="ML")
ajusteSqrt## Series: train2.ts
## Regression with ARIMA(1,0,0)(1,0,1)[12] errors
##
## Coefficients:
## ar1 sar1 sma1 intercept semanaSanta dt bisiesto
## 0.4832 0.9583 -0.6101 117.2187 1.9221 -0.5387 0.8696
## s.e. 0.0915 0.0287 0.1206 7.9647 0.3536 0.1298 3.2162
## festivo price
## 4.4400 -3.3974
## s.e. 1.7961 1.2340
##
## sigma^2 estimated as 27.68: log likelihood=-298.3
## AIC=616.59 AICc=619.18 BIC=642.24summary(ajusteSqrt)## Series: train2.ts
## Regression with ARIMA(1,0,0)(1,0,1)[12] errors
##
## Coefficients:
## ar1 sar1 sma1 intercept semanaSanta dt bisiesto
## 0.4832 0.9583 -0.6101 117.2187 1.9221 -0.5387 0.8696
## s.e. 0.0915 0.0287 0.1206 7.9647 0.3536 0.1298 3.2162
## festivo price
## 4.4400 -3.3974
## s.e. 1.7961 1.2340
##
## sigma^2 estimated as 27.68: log likelihood=-298.3
## AIC=616.59 AICc=619.18 BIC=642.24
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.4333323 5.008137 3.950826 -0.6640638 3.855063 0.6270734
## ACF1
## Training set 0.008130705
Vemos en este caso que también conseguimos capturar todas las barras de correlaciones dentro de las bandas de Barlett, mejorando con los regresores externos el modelo expuesto anteriormente.
Tercer Modelo
Queremos probar también el otro modelo Arima sin transformaciones que habíamos planteado antes (0,1,1)x(0,1,1) con las variables exógenas. Esa era la intención, pero nos encontramos con un error persistente de R en el que nos dice “non-finite value supplied by optim”. Curiosamente este error sólo desaparece cuando quitamos la diferencia estacional: (1, 1, 1)x(1, 0, 1). Sin embargo esa diferencia estacional, una vez hemos decidido añadir SAr1, es necesaria… Por eso nos quedamos de momento con ese modelo tal cual, sin incluir las variables exógenas. Por otro lado, en el Jupyter Notebook adjunto en la carpeta que contiene este archivo vamos a intentar resolver este problema que nos hemos encontrado con Python, tratando de ver si éste puede optimizar la función que le especificamos.
En primer lugar, debe decirse que Python también tuvo problemas para optimizar el modelo SARIMAX que proponíamos, pero está algo mejor preparado el paquete statsmodels y tiene un método que le permite cambiar la fórmula utilizada para optimizar, sin utilizar gradientes ni hessianos, lo cual estaba causando problemas por algún motivo. Sin embargo, el ajuste es bastante malo (ver las estadísticas del Jupyter Notebook, en ellas se puede apreciar que casi ningún coeficiente sale significativo), una muestra de que existe un problema matemático al juntar el modelo con esa diferencia estacional con los regresores externos; al no poder utilizar los mejores métodos de optimización que son los que vienen por defecto, la estimación de los parámetros estará realizada de una forma más defectuosa, y por lo tanto el modelo no es tan bueno. Una muestra de esto es que, como vemos en los gráficos superiores, existe una autocorrelación parcial muy significativa con un período muy anterior, de muchos años… Esto es raro y es un comportamiento de los residuos que no hemos observado en el resto de modelos.
Cuarto Modelo.
Probaremos ahora con el modelo de \(y^\lambda\). De nuevo, quitamos la diferencia estacional para poder realizar los cálculos, pues con ella R no es capaz de optimizar la función y por lo tanto no puede estimar los parámetros.
ajusteLambda= Arima(train3.ts,
order=c(1,1,1),
seasonal=list(order=c(1,0,1), period=12),
xreg=xregs,
method="ML")
ajusteLambda## Series: train3.ts
## Regression with ARIMA(1,1,1)(1,0,1)[12] errors
##
## Coefficients:
## ar1 ma1 sar1 sma1 semanaSanta dt bisiesto
## 0.3718 -0.9443 0.9765 -0.6683 0.2617 -0.0714 0.1229
## s.e. 0.1049 0.0311 0.0184 0.1129 0.0476 0.0172 0.4205
## festivo price
## 0.5842 0.1647
## s.e. 0.2670 0.2538
##
## sigma^2 estimated as 0.4488: log likelihood=-101.77
## AIC=223.55 AICc=226.17 BIC=249.09summary(ajusteLambda)## Series: train3.ts
## Regression with ARIMA(1,1,1)(1,0,1)[12] errors
##
## Coefficients:
## ar1 ma1 sar1 sma1 semanaSanta dt bisiesto
## 0.3718 -0.9443 0.9765 -0.6683 0.2617 -0.0714 0.1229
## s.e. 0.1049 0.0311 0.0184 0.1129 0.0476 0.0172 0.4205
## festivo price
## 0.5842 0.1647
## s.e. 0.2670 0.2538
##
## sigma^2 estimated as 0.4488: log likelihood=-101.77
## AIC=223.55 AICc=226.17 BIC=249.09
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.09179785 0.6340547 0.5019316 -0.5206396 2.413343 0.5958316
## ACF1
## Training set 0.003080306
En este caso los residuos sí presentan aún correlaciones significativas estadísticamente, por lo que no hemos conseguido mejorar demasiado en este caso introduciendo variables exógenas. Como vemos que hay algunas de estas variables que no son significativas (como el precio, por ejemplo), en este caso vamos a darle sólo algunas columnas de la matriz de regresores externos, las que pertenecen a las variables significativas.
## Series: train3.ts
## Regression with ARIMA(1,1,1)(1,0,1)[12] errors
##
## Coefficients:
## ar1 ma1 sar1 sma1 semanaSanta dt festivo
## 0.3785 -0.9489 0.9744 -0.6598 0.2611 -0.0716 0.5841
## s.e. 0.1052 0.0312 0.0194 0.1129 0.0476 0.0173 0.2627
##
## sigma^2 estimated as 0.4439: log likelihood=-102.04
## AIC=220.07 AICc=221.75 BIC=240.5## Series: train3.ts
## Regression with ARIMA(1,1,1)(1,0,1)[12] errors
##
## Coefficients:
## ar1 ma1 sar1 sma1 semanaSanta dt festivo
## 0.3785 -0.9489 0.9744 -0.6598 0.2611 -0.0716 0.5841
## s.e. 0.1052 0.0312 0.0194 0.1129 0.0476 0.0173 0.2627
##
## sigma^2 estimated as 0.4439: log likelihood=-102.04
## AIC=220.07 AICc=221.75 BIC=240.5
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.08939405 0.6378955 0.5088941 -0.5121413 2.445117 0.6040966
## ACF1
## Training set 0.002436342Hay algo que venimos observando en la mayoría de los modelos, y es que cuando no introducimos una diferencia estacional el modelo lo pide a gritos, pues \(SAR1 \pm 2 \cdot sd(SAR1) > 1\).
Una vez hemos quitado las variables exógenas que estaban introduciendo ruido parece que el modelo se ajusta mucho mejor, capturando casi todas las autocorrelaciones tanto simples como parciales dentro de las barras de Barlett.
Quinto Modelo: Ajuste Automático
Probaremos también el ajuste automático del paquete forecast, para compararlo con los modelos que acabamos de construir.
## Series: train.ts
## Regression with ARIMA(1,0,0)(1,0,0)[12] errors
##
## Coefficients:
## ar1 sar1 intercept semanaSanta dt bisiesto
## 0.4537 0.6961 12449.662 378.4617 -111.4025 -287.0549
## s.e. 0.0927 0.0775 1882.091 72.3948 34.9729 702.3373
## festivo price
## 840.6330 -467.3436
## s.e. 306.3595 317.3659
##
## sigma^2 estimated as 1551111: log likelihood=-820.35
## AIC=1658.7 AICc=1660.79 BIC=1681.77
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -58.6853 1192.414 918.0954 -1.832455 8.756608 0.7034693
## ACF1
## Training set 0.03301758##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 4.5372e-01 9.2732e-02 4.8928 9.942e-07 ***
## sar1 6.9607e-01 7.7462e-02 8.9860 < 2.2e-16 ***
## intercept 1.2450e+04 1.8821e+03 6.6148 3.720e-11 ***
## semanaSanta 3.7846e+02 7.2395e+01 5.2277 1.716e-07 ***
## dt -1.1140e+02 3.4973e+01 -3.1854 0.001446 **
## bisiesto -2.8705e+02 7.0234e+02 -0.4087 0.682750
## festivo 8.4063e+02 3.0636e+02 2.7439 0.006071 **
## price -4.6734e+02 3.1737e+02 -1.4726 0.140867
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## ar1 sar1 intercept semanaSanta dt
## ar1 1.000000000 0.02869996 -0.17252556 0.034373306 0.03105312
## sar1 0.028699963 1.00000000 -0.22470453 0.015800651 0.06786113
## intercept -0.172525564 -0.22470453 1.00000000 -0.032907182 -0.02962391
## semanaSanta 0.034373306 0.01580065 -0.03290718 1.000000000 0.16617791
## dt 0.031053119 0.06786113 -0.02962391 0.166177908 1.00000000
## bisiesto -0.032525076 -0.13584247 0.05270572 0.038990950 0.06340133
## festivo -0.002152932 -0.08497891 -0.11188607 0.048950588 -0.03295701
## price 0.178976587 0.23922635 -0.93846931 0.008815508 0.02856165
## bisiesto festivo price
## ar1 -0.03252508 -0.002152932 0.178976587
## sar1 -0.13584247 -0.084978910 0.239226352
## intercept 0.05270572 -0.111886071 -0.938469314
## semanaSanta 0.03899095 0.048950588 0.008815508
## dt 0.06340133 -0.032957006 0.028561651
## bisiesto 1.00000000 0.105985726 -0.083067380
## festivo 0.10598573 1.000000000 -0.043296371
## price -0.08306738 -0.043296371 1.000000000## Retard p-value
## [1,] 6 0.9712270
## [2,] 12 0.5183626
## [3,] 18 0.6650078
## [4,] 24 0.6116338
## [5,] 30 0.5438551
## [6,] 36 0.5760761
## [7,] 42 0.6409670
## [8,] 48 0.6901314
En este caso el precio parece no ser significativo, tendríamos que coger \(\alpha = 0.15\) para que éste fuera significativo. Esto podría estar introduciendo algo de ruido… El Stepwise lo coge porque seguramente el AIC baje al introducirlo, pero no nos convence.
ajusteAutomatico2<-auto.arima(train.ts,max.d=1, max.D=1, max.p=1, max.q=1, max.P=1, max.Q=1,
seasonal=TRUE,ic="aic",xreg=xregs[ , 1:4],stepwise=TRUE)
summary(ajusteAutomatico2)## Series: train.ts
## Regression with ARIMA(1,1,1)(1,0,0)[12] errors
##
## Coefficients:
## ar1 ma1 sar1 semanaSanta dt bisiesto festivo
## 0.4309 -0.9603 0.7183 378.2619 -109.6057 -306.8391 833.8064
## s.e. 0.0998 0.0254 0.0748 71.3012 34.5435 688.3324 319.6909
##
## sigma^2 estimated as 1534906: log likelihood=-812.37
## AIC=1640.74 AICc=1642.41 BIC=1661.17
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -155.7495 1186.169 917.4091 -2.598247 8.826862 0.7029434
## ACF1
## Training set 0.01743244coeftest(ajusteAutomatico2)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.430888 0.099837 4.3159 1.589e-05 ***
## ma1 -0.960282 0.025393 -37.8173 < 2.2e-16 ***
## sar1 0.718279 0.074777 9.6056 < 2.2e-16 ***
## semanaSanta 378.261909 71.301206 5.3051 1.126e-07 ***
## dt -109.605727 34.543479 -3.1730 0.001509 **
## bisiesto -306.839073 688.332416 -0.4458 0.655762
## festivo 833.806369 319.690888 2.6082 0.009103 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cor.arma(ajusteAutomatico2)## ar1 ma1 sar1 semanaSanta dt
## ar1 1.000000000 -0.32387322 -0.007077443 0.03629576 0.02554625
## ma1 -0.323873218 1.00000000 -0.202175178 -0.01596965 -0.01307942
## sar1 -0.007077443 -0.20217518 1.000000000 0.02190270 0.06614748
## semanaSanta 0.036295761 -0.01596965 0.021902700 1.00000000 0.16615284
## dt 0.025546254 -0.01307942 0.066147478 0.16615284 1.00000000
## bisiesto -0.031243024 0.04894940 -0.115393743 0.03544888 0.06710558
## festivo -0.003781187 0.02558829 -0.069071152 0.04473960 -0.03011573
## bisiesto festivo
## ar1 -0.03124302 -0.003781187
## ma1 0.04894940 0.025588287
## sar1 -0.11539374 -0.069071152
## semanaSanta 0.03544888 0.044739603
## dt 0.06710558 -0.030115726
## bisiesto 1.00000000 0.097608925
## festivo 0.09760893 1.000000000Box.test.2(residuals(ajusteAutomatico2),
nlag = c(6,12,18,24,30,36,42,48),
type="Ljung-Box")## Retard p-value
## [1,] 6 0.9244382
## [2,] 12 0.2802904
## [3,] 18 0.4741342
## [4,] 24 0.4239821
## [5,] 30 0.3697063
## [6,] 36 0.4117417
## [7,] 42 0.4889451
## [8,] 48 0.5470928acf(ajusteAutomatico2$residuals, lag.max = 48, xlab = "Retardo",
main= "Función de autocorrelación simple")
pacf(ajusteAutomatico2$residuals, lag.max = 48, xlab = "Retardo",
main = "Función de autocorrelación parcial")
Vemos que no cambia demasiado con respecto a cuando metemos price como regresor externo. Nos quedamos con el primero, pues tiene un MAPE algo menor.
Identificación de Outliers
En esta sección nos vamos a encargar de detectar los outliers y las anomalías de los modelos que hemos empleado en R, pues el que hemos empleado en Python tendría que ser tratado con paquetes de Python y con otro tipo de procedimientos. Nos centramos por tanto en los siguientes modelos:
Primer Modelo: SARIMAX(0, 0, 1)x(1, 0, 1)[12]
Segundo Modelo: \(\sqrt{y}\) SARIMAX(1, 0, 0)x(1, 0, 1)[12]
Tercer Modelo: \(y^\lambda\) SARIMAX (0, 1, 1)x(1, 0, 1)[12]
Cuarto Modelo: Ajuste Automático SARIMAX(1,1,1)x(1,0,0)[12]
En este apartado utilizaremos siglas para referirnos a los diferentes tipos de anomalías o outliers que existen en una serie temporal, para abreviar el texto. Por ello, aquí vamos a especificar de qué tipo de outliers hablamos cuando utilizamos cada una de estas siglas.
AO: Additive Outliers
IO: Innovative Outliers: estos son un tipo de outliers que dependen del modelo que se esté utilizando, en contraste con el resto. Más información sobre este tema, así sobre el resto de tipos de outliers aquí.
LS: Level Shift.
TC: Transitory Change.
El propio nombre de cada tipo de outlier es suficientemente informativo, ya que nos dice cuál es el tipo de cambio que sufre la serie en esa observación que consideramos outlier.
Identificación de los Outliers del Primer Modelo
## Warning: package 'tsoutliers' was built under R version 3.5.2## ind type coefhat tstat abststat date
## 1 43 TC 3221.6158 3.998365 3.998365 2004-07-01
## 2 67 LS -992.9195 -3.040895 3.040895 2006-07-01
## 3 70 LS -1007.7600 -3.036758 3.036758 2006-10-01
## 4 81 AO -2837.9402 -3.356161 3.356161 2007-09-01Vaya, parece que si encontramos algunos patrones de outliers. En concreto, vemos que tenemos un outlier repetido, aunque de diferente tipo (el primero es de tipo TC y el segundo de tipo LS) el 1 de Julio de 2004 y 2006. Por otro lado, el 1 de Octubre de 2006 tenemos otro del tipo LS y el 1 de Septiembre del 2007 tenemos uno del tipo AO.
library(ggplot2)
v = rep(0, nrow(train))
v[c(40:46, 66:68, 69:71, 78:84)] <- 1
inds = diff(c(0, v))
start = train$date[inds == 1]
end = train$date[inds==-1]
if(length(start) > length(end)) end <- c(end, tail(train$date, 1))
rects = data.frame(start=start, end = end, group=seq_along(start))
ggplot(data = train, aes(date, users)) +
theme_minimal() +
geom_line(lty=2, color="steelblue", lwd=1.1) +
geom_point() +
geom_rect(data = rects, inherit.aes=FALSE, aes(xmin=start, xmax=end, ymin=min(train$users),
ymax=max(train$users), group=group), color="transparent", fill="orange", alpha=0.3) +
geom_vline(aes(xintercept=as.Date("2007-09-01"))) +
geom_vline(aes(xintercept=as.Date("2004-07-01"))) +
geom_vline(aes(xintercept=as.Date("2006-07-01"))) +
geom_vline(aes(xintercept=as.Date("2006-10-01"))) +
ggtitle("Outliers of the Train Set with the First Model")
En el gráfico superior podemos ver los diferentes outliers que presenta el primer modelo.
## ind type coefhat tstat abststat date
## 1 43 TC 3221.6158 3.998365 3.998365 2004-07-01
## 2 81 AO -2837.9402 -3.356161 3.356161 2007-09-01
## 3 67 LS -992.9195 -3.040895 3.040895 2006-07-01
## 4 70 LS -1007.7600 -3.036758 3.036758 2006-10-01Ahora ajustaremos de nuevo el modelo teniendo en cuenta los outliers. Vemos que las variables LS67 y LS70 no son estadísticamente significativas, por lo que las eliminamos. Eliminamos también bisiesto que no es significativa.
ajuste5def = Arima(train.ts,
order = c(0,0,1),
seasonal = list(order = c(1,0,1), period = 12),
xreg = RegressorsMasOutliers[, c(1, 2, 4,5,6,7)],
method = "ML")
ajuste5def## Series: train.ts
## Regression with ARIMA(0,0,1)(1,0,1)[12] errors
##
## Coefficients:
## ma1 sar1 sma1 intercept semanaSanta dt festivo
## 0.4718 0.9553 -0.5802 13902.622 356.4204 -127.0938 862.4788
## s.e. 0.0929 0.0296 0.1263 1197.981 61.6411 26.0868 335.2374
## price TC43 AO81
## -737.9075 3133.9434 -2871.3002
## s.e. 186.0039 757.3415 786.1132
##
## sigma^2 estimated as 951818: log likelihood=-799.23
## AIC=1620.45 AICc=1623.6 BIC=1648.66summary(ajuste5def)## Series: train.ts
## Regression with ARIMA(0,0,1)(1,0,1)[12] errors
##
## Coefficients:
## ma1 sar1 sma1 intercept semanaSanta dt festivo
## 0.4718 0.9553 -0.5802 13902.622 356.4204 -127.0938 862.4788
## s.e. 0.0929 0.0296 0.1263 1197.981 61.6411 26.0868 335.2374
## price TC43 AO81
## -737.9075 3133.9434 -2871.3002
## s.e. 186.0039 757.3415 786.1132
##
## sigma^2 estimated as 951818: log likelihood=-799.23
## AIC=1620.45 AICc=1623.6 BIC=1648.66
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -94.48002 923.4014 726.6795 -1.659199 6.872925 0.5568013
## ACF1
## Training set 0.02273563Identificación de Outliers del Segundo Modelo
## ind type coefhat tstat abststat date
## 1 81 AO -14.58462 -3.415948 3.415948 2007-09-01En este caso solamente tenemos un outlier del tipo AO el 1 de Septiembre de 2007.
En este caso eliminamos la variable bisiesto pues no es significativa, y ajustamos el siguiente modelo:
## Series: train2.ts
## Regression with ARIMA(1,0,0)(1,0,1)[12] errors
##
## Coefficients:
## ar1 sar1 sma1 intercept semanaSanta dt festivo
## 0.5230 0.9618 -0.5908 116.1457 1.8731 -0.6526 4.2460
## s.e. 0.0901 0.0260 0.1241 8.2292 0.3164 0.1201 1.7845
## price AO81
## -3.1451 -15.4789
## s.e. 1.2509 3.9725
##
## sigma^2 estimated as 23.6: log likelihood=-291.45
## AIC=602.89 AICc=605.48 BIC=628.54## Series: train2.ts
## Regression with ARIMA(1,0,0)(1,0,1)[12] errors
##
## Coefficients:
## ar1 sar1 sma1 intercept semanaSanta dt festivo
## 0.5230 0.9618 -0.5908 116.1457 1.8731 -0.6526 4.2460
## s.e. 0.0901 0.0260 0.1241 8.2292 0.3164 0.1201 1.7845
## price AO81
## -3.1451 -15.4789
## s.e. 1.2509 3.9725
##
## sigma^2 estimated as 23.6: log likelihood=-291.45
## AIC=602.89 AICc=605.48 BIC=628.54
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.4389988 4.625089 3.669041 -0.6284198 3.555108 0.5823486
## ACF1
## Training set -0.004087827##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.523007 0.090058 5.8075 6.343e-09 ***
## sar1 0.961818 0.026033 36.9464 < 2.2e-16 ***
## sma1 -0.590811 0.124114 -4.7602 1.934e-06 ***
## intercept 116.145686 8.229197 14.1139 < 2.2e-16 ***
## semanaSanta 1.873100 0.316407 5.9199 3.221e-09 ***
## dt -0.652622 0.120132 -5.4325 5.556e-08 ***
## festivo 4.245998 1.784493 2.3794 0.01734 *
## price -3.145081 1.250873 -2.5143 0.01193 *
## AO81 -15.478855 3.972484 -3.8965 9.759e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Identificación de Outliers del Tercer Modelo
## ind type coefhat tstat abststat date
## 1 81 AO -1.924749 -3.610469 3.610469 2007-09-01Sale exactamente el mismo outlier que en el caso anterior (también había aparecido en el primero). Esto contradice lo que se mencionaba en el link que se encuentra más arriba, en el que se decía que los outliers del tipo IO no dependen del modelo a diferencia de los demás. Por lo que estamos observando, el outlier que no depende del modelo es el AO de 2007-09-01.
## $title
## [1] "Outliers of the Train Set with the Third Model"
##
## attr(,"class")
## [1] "labels"En este caso eliminamos bisiesto y price, pues no son significativas.
ajusteLambdaDef = Arima(train3.ts,
order = c(0,1,1),
seasonal = list(order = c(1,0,1), period = 12),
xreg = RegressorsMasOutliers3[, c(1, 2, 4, 6)],
method = "ML")
ajusteLambdaDef## Series: train3.ts
## Regression with ARIMA(0,1,1)(1,0,1)[12] errors
##
## Coefficients:
## ma1 sar1 sma1 semanaSanta dt festivo AO81
## -0.5740 0.9820 -0.6699 0.2548 -0.0900 0.5410 -2.2398
## s.e. 0.1552 0.0147 0.1198 0.0478 0.0178 0.2848 0.5493
##
## sigma^2 estimated as 0.4275: log likelihood=-101.68
## AIC=219.36 AICc=221.03 BIC=239.79summary(ajusteLambdaDef)## Series: train3.ts
## Regression with ARIMA(0,1,1)(1,0,1)[12] errors
##
## Coefficients:
## ma1 sar1 sma1 semanaSanta dt festivo AO81
## -0.5740 0.9820 -0.6699 0.2548 -0.0900 0.5410 -2.2398
## s.e. 0.1552 0.0147 0.1198 0.0478 0.0178 0.2848 0.5493
##
## sigma^2 estimated as 0.4275: log likelihood=-101.68
## AIC=219.36 AICc=221.03 BIC=239.79
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.02216439 0.6260056 0.508905 -0.1635231 2.43026 0.6041096
## ACF1
## Training set 0.1444174coeftest(ajusteLambdaDef)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 -0.574004 0.155191 -3.6987 0.0002167 ***
## sar1 0.981959 0.014662 66.9715 < 2.2e-16 ***
## sma1 -0.669914 0.119799 -5.5920 2.245e-08 ***
## semanaSanta 0.254823 0.047787 5.3324 9.691e-08 ***
## dt -0.089971 0.017800 -5.0545 4.315e-07 ***
## festivo 0.541018 0.284765 1.8999 0.0574490 .
## AO81 -2.239846 0.549277 -4.0778 4.546e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Identificación de Outliers del Cuarto Modelo
## ind type coefhat tstat abststat date
## 1 26 AO 2395.624 3.290378 3.290378 2003-02-01
## 2 40 TC 2295.855 3.029250 3.029250 2004-04-01
## 3 43 TC 2959.651 3.905092 3.905092 2004-07-01
## 4 81 AO -2919.068 -4.009327 4.009327 2007-09-01En este caso observamos unos outliers algo distintos, en concreto el del 1 de Frebrero del 2003.
## Series: train.ts
## Regression with ARIMA(0,0,1)(1,0,0)[12] errors
##
## Coefficients:
## ma1 sar1 intercept semanaSanta dt bisiesto
## 0.5725 0.7442 12839.046 323.7978 -144.1346 985.1539
## s.e. 0.1154 0.0742 1398.947 50.7277 29.1279 564.1778
## festivo price AO26 TC43 AO81
## 862.9701 -555.2582 2915.7790 3174.1653 -3092.7717
## s.e. 257.3901 234.4697 710.8307 718.1577 664.0621
##
## sigma^2 estimated as 1047564: log likelihood=-800.79
## AIC=1625.58 AICc=1629.34 BIC=1656.35
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -67.8566 963.0838 778.422 -1.539647 7.378763 0.5964477
## ACF1
## Training set 0.000512715##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 5.7247e-01 1.1541e-01 4.9602 7.041e-07 ***
## sar1 7.4415e-01 7.4153e-02 10.0354 < 2.2e-16 ***
## intercept 1.2839e+04 1.3989e+03 9.1777 < 2.2e-16 ***
## semanaSanta 3.2380e+02 5.0728e+01 6.3831 1.736e-10 ***
## dt -1.4413e+02 2.9128e+01 -4.9483 7.485e-07 ***
## bisiesto 9.8515e+02 5.6418e+02 1.7462 0.0807803 .
## festivo 8.6297e+02 2.5739e+02 3.3528 0.0008001 ***
## price -5.5526e+02 2.3447e+02 -2.3681 0.0178775 *
## AO26 2.9158e+03 7.1083e+02 4.1019 4.097e-05 ***
## TC43 3.1742e+03 7.1816e+02 4.4199 9.876e-06 ***
## AO81 -3.0928e+03 6.6406e+02 -4.6574 3.203e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## ma1 sar1 intercept semanaSanta
## ma1 1.000000000 -0.20437116 -0.135490559 -0.203795824
## sar1 -0.204371158 1.00000000 -0.206725827 0.068018944
## intercept -0.135490559 -0.20672583 1.000000000 -0.013153707
## semanaSanta -0.203795824 0.06801894 -0.013153707 1.000000000
## dt -0.005670653 0.06776805 -0.045267277 0.203753673
## bisiesto 0.261945821 -0.15536832 0.048395787 -0.091535287
## festivo -0.018197633 -0.09188222 -0.124342996 0.032923276
## price 0.146219717 0.22662111 -0.928634546 -0.007164621
## AO26 0.095239154 -0.09570366 0.123411439 -0.174263687
## TC43 0.220047765 -0.12254738 -0.001739071 0.006132956
## AO81 -0.068885427 0.04353456 -0.054292399 0.077101984
## dt bisiesto festivo price
## ma1 -0.005670653 0.26194582 -0.018197633 0.146219717
## sar1 0.067768054 -0.15536832 -0.091882220 0.226621113
## intercept -0.045267277 0.04839579 -0.124342996 -0.928634546
## semanaSanta 0.203753673 -0.09153529 0.032923276 -0.007164621
## dt 1.000000000 -0.08464818 -0.067611369 0.052000771
## bisiesto -0.084648181 1.00000000 0.111077678 -0.083761033
## festivo -0.067611369 0.11107768 1.000000000 -0.050783867
## price 0.052000771 -0.08376103 -0.050783867 1.000000000
## AO26 -0.146432322 0.40544776 0.070151245 -0.148718697
## TC43 0.060416959 0.09769199 0.006446041 -0.010730443
## AO81 0.277853894 -0.02177463 -0.002048460 0.051944791
## AO26 TC43 AO81
## ma1 0.095239154 0.220047765 -0.068885427
## sar1 -0.095703658 -0.122547377 0.043534558
## intercept 0.123411439 -0.001739071 -0.054292399
## semanaSanta -0.174263687 0.006132956 0.077101984
## dt -0.146432322 0.060416959 0.277853894
## bisiesto 0.405447761 0.097691991 -0.021774634
## festivo 0.070151245 0.006446041 -0.002048460
## price -0.148718697 -0.010730443 0.051944791
## AO26 1.000000000 0.009820674 -0.045422480
## TC43 0.009820674 1.000000000 0.001388494
## AO81 -0.045422480 0.001388494 1.000000000## Retard p-value
## [1,] 6 0.10212090
## [2,] 12 0.03233052
## [3,] 18 0.01566603
## [4,] 24 0.00889139
## [5,] 30 0.00779644
## [6,] 36 0.01289081
## [7,] 42 0.01397820
## [8,] 48 0.02218604
Vemos que la serie sigue presentando anomalías. Además, el ajuste automático ha cambiado, y ahora estima un SARIMAX(0, 0, 1)x(1, 0, 0).
Análisis de la Varianza
Primer Modelo
La pregunta que nos hacemos antes de ponernos a predecir es la siguiente, ¿admitiría alguno de los modelos que hemos visto arriba el ajuste de un modelo GARCH para modelizar la varianza de la serie?.
## Retard p-value
## [1,] 6 0.6938020
## [2,] 12 0.8881881
## [3,] 18 0.7590623
## [4,] 24 0.8721452
## [5,] 30 0.9708096
## [6,] 36 0.9898367
## [7,] 42 0.9969951
## [8,] 48 0.9829229En este caso los residuos presentan ruido blanco, por lo tanto no vamos a poder utilizar un modelo GARCH.
Segundo Modelo
## Retard p-value
## [1,] 6 0.05882292
## [2,] 12 0.12071514
## [3,] 18 0.02685252
## [4,] 24 0.07757110
## [5,] 30 0.07584400
## [6,] 36 0.17743292
## [7,] 42 0.21977585
## [8,] 48 0.28055270En este caso el ruido blanco no está tan claro, por lo que sí podría tener sentido utilizar un modelo GARCH para modelizar la varianza de la serie. Sin embargo, la modelización de la varianza de la serie temporal va a ser llevada acabo únicamente con la serie ajustada de manera automática, pues lo vamos a utilizar para contrastar un “Machine Learning Approach” versus los ajustes manuales que hace el humano con su conocimiento de estadística y su detección de patrones (al final lo que hacemos al analizar los gráficos ACF y PACF y ajustando de nuevo los parámetros es detectar patrones en la serie temporal y jugar con ellos).
Tercer Modelo
## Retard p-value
## [1,] 6 0.02051878
## [2,] 12 0.09414880
## [3,] 18 0.26736440
## [4,] 24 0.30652493
## [5,] 30 0.21065611
## [6,] 36 0.31032755
## [7,] 42 0.52150333
## [8,] 48 0.58652776Algo menos que en el anterior caso, quizás sí se podría utilizar un GARCH.
Cuarto Modelo
## Retard p-value
## [1,] 6 0.3452732
## [2,] 12 0.5060421
## [3,] 18 0.3245441
## [4,] 24 0.4997114
## [5,] 30 0.7085126
## [6,] 36 0.8029283
## [7,] 42 0.2567279
## [8,] 48 0.3872969En este caso queda claro que el proceso es de ruido blanco y por lo tanto no parece tener mucho sentido ajustar un modelo GARCH. Lo que vamos a hacer en la siguiente sección será ajustar primero los 3 modelos manuales únicamente para el modelo de la media, y por otro lado, con una función diferente, vamos a realizar un ajuste automático en el que simultáneamente predigamos el valor medio, así como los valores superiores e inferiores, obtenidos estos mediante la predicción de la desviación estándar con un modelo GARCH(1,1).
PREDICCIONES, EVALUACIÓN DE MODELOS Y RESULTADOS
Time Series Approach
get_predictions = function(d = datos, nendtrain=96, var='users', datecol='date', lam=0.6566137){
########################################################################################################################
# This function will calculate the predictions, in periods of 12 months, for the cinema users in Spain #
# with the different models presented. #
# Parameters: #
# - d: data #
# - nendtrain: index of the last observation for the pure train part, before starting to predict (moving window)#
# - var: the variable of interest inside the dataset for univariate time series #
# - datecol: the name or index of the date variable in the dataset #
# - lam: the lambda value for Box-Cox transformation as returned by car::powerTransform() #
# Returns: #
# - preds_df: data frame with predictions and real value. #
#######################################################################################################################
preds1 = c()
preds2 = c()
preds3 = c()
for (i in seq(nendtrain, nrow(d) - 12, 12)){
#dates = d[1:i, datecol]
###########################
# TRAIN & TEST #
# SPLITS #
###########################
tr1 = data.frame(d[1:i,])
tr2 = data.frame(d[1:i, ])
tr2[, var] = sqrt(tr2[, var])
tr3 = data.frame(d[1:i, ])
tr3[, var] = tr3[,var]**lam
te1 = data.frame(d[(i+1):(i+12), ])
te2 = data.frame(d[(i+1):(i+12), ])
te2[,var] = sqrt(te2[,var])
te3 = data.frame(d[(i+1):(i+12), ])
te3[,var] = te3[,var]**lam
#####################
# COMPUTING #
# EXTERNAL VARIABLES#
#####################
xregs = as.matrix(calculoExplicativasCalendario(tr1[,datecol],0)[ , c("semanaSanta", "dt", "bisiesto", "festivo", "price")])
newxregs = as.matrix(calculoExplicativasCalendario(te1$date, 0)[,c("semanaSanta","dt","bisiesto","festivo","price")])
tr1.ts = ts(tr1[,var], start=c(year(tr1[1, datecol]), month(tr1[1, datecol])),
end=c(year(tr1[nrow(tr1), datecol]), month(tr1[nrow(tr1), datecol])),
frequency=12)
tr2.ts = ts(tr2[,var], start=c(year(tr2[1, datecol]), month(tr2[1, datecol])),
end=c(year(tr2[nrow(tr2), datecol]), month(tr2[nrow(tr2), datecol])),
frequency=12)
tr3.ts = ts(tr3[,var], start=c(year(tr3[1, datecol]), month(tr3[1, datecol])),
end=c(year(tr3[nrow(tr3), datecol]), month(tr3[nrow(tr3), datecol])),
frequency=12)
te1.ts = ts(te1[, var], start=c(year(te1[1, datecol]), month(te1[1,datecol])),
end = c(year(te1[nrow(te1), datecol]), month(te1[nrow(te1),datecol])),
frequency=12)
te2.ts = ts(te2[, var], start=c(year(te2[1, datecol]), month(te2[1,datecol])),
end = c(year(te2[nrow(te2), datecol]), month(te2[nrow(te2),datecol])),
frequency=12)
te3.ts = ts(te3[, var], start=c(year(te3[1, datecol]), month(te3[1,datecol])),
end = c(year(te3[nrow(te3), datecol]), month(te3[nrow(te3),datecol])),
frequency=12)
###############################################################################################################
# FITTING, OUTLIERS CALCULATION, EXTERNAL VARIABLES RECALCULATED WITH OUTLIERS; PREDICTION FOR EACH MODEL #
###############################################################################################################
#################
# 1ST MODEL #
#################
ajuste1 <- Arima(tr1.ts,
order = c(0,0,1),
seasonal = list(order = c(1,0,1), period = 12),
xreg = xregs,
method = "ML")
listaOutliers1 <- locate.outliers(ajuste1$residuals,
pars = coefs2poly(ajuste1),
types = c("AO", "LS", "TC"),cval=3)
outliers1 <- outliers(c(as.character(listaOutliers1$type)), c(listaOutliers1$ind))
outliersVariables1 <- outliers.effects(outliers1, length(ajuste1$residuals))
RegressorsMasOutliers1 <- cbind(xregs,outliersVariables1)
ajuste1def=Arima(tr1.ts,
order = c(0,0,1),
seasonal = list(order = c(1,0,1), period = 12),
xreg = RegressorsMasOutliers1[, -3],
method = "ML")
outliersNew1 <- outliersVariables1[(i-11):i,]
newxregs1 <- as.matrix(cbind(newxregs,outliersNew1))
prediccion1 <- predict(ajuste1def, n.ahead=12,
newxreg=newxregs1[, -3])
preds1 = c(preds1, prediccion1$pred)
#################
# 2ND MODEL #
#################
ajuste2 = Arima(tr2.ts,
order=c(1,0,0),
seasonal=list(order=c(1,0,1), period=12),
xreg=xregs,
method="ML")
listaOutliers2 <- locate.outliers(ajuste2$residuals,
pars = coefs2poly(ajuste2),
types = c("AO", "LS", "TC"),cval=3)
outliers2 <- outliers(c(as.character(listaOutliers2$type)), c(listaOutliers2$ind))
tryCatch(
{outliersVariables2 <- outliers.effects(outliers2, length(ajuste2$residuals))
RegressorsMasOutliers2 <- cbind(xregs,outliersVariables2)
},
error=function(e){RegressorsMasOutliers2 = as.matrix(xregs)}
)
if(nrow(RegressorsMasOutliers2 != length(tr2.ts))){
print("check this problem again")
a=0
RegressorsMasOutliers2=as.matrix(xregs)
}else{a = 1}
ajuste2def = Arima(tr2.ts,
order = c(1,0,0),
seasonal = list(order = c(1,0,1), period = 12),
xreg = RegressorsMasOutliers2[, -3],
method = "ML")
if(a == 0){
newxregs2 <- as.matrix(newxregs[,-3])
prediccion2 <- predict(ajuste2def, n.ahead=12,
newxreg=newxregs2)
preds2 = c(preds2, prediccion2$pred)
}else{
outliersNew2 <- outliersVariables2[(i-11):i,]
newxregs2 <- as.matrix(cbind(newxregs,outliersNew2))
prediccion2 <- predict(ajuste2def, n.ahead=12,
newxreg=newxregs2[, -3])
preds2 = c(preds2, prediccion2$pred)
}
###############
# 3RD MODEL #
###############
ajuste3 = Arima(tr3.ts,
order=c(1,1,1),
seasonal=list(order=c(1,0,1), period=12),
xreg=xregs[, -c(3,5)],
method="ML")
listaOutliers3 <- locate.outliers(ajuste3$residuals,
pars = coefs2poly(ajuste3),
types = c("AO", "LS", "TC"),cval=3)
outliers3 <- outliers(c(as.character(listaOutliers3$type)), c(listaOutliers3$ind))
outliersVariables3 <- outliers.effects(outliers3, length(ajuste3$residuals))
RegressorsMasOutliers3 <- cbind(xregs,outliersVariables3)
ajuste3def=Arima(tr3.ts,
order = c(0,1,1),
seasonal = list(order = c(1,0,1), period = 12),
xreg = RegressorsMasOutliers3[, -c(3,5)],
method = "ML")
outliersNew3 <- outliersVariables3[(i-11):i,]
newxregs3 <- as.matrix(cbind(newxregs,outliersNew3))
prediccion3 <- predict(ajuste3def, n.ahead=12,
newxreg=newxregs3[, -c(3,5)])
preds3 = c(preds3, prediccion3$pred)
}
if(length(preds1)!=length(preds2) | length(preds1) != length(preds3)){
return(list(preds1, preds2, preds3))
}
real = data.frame(d)[(nendtrain+1):nrow(d), var]
preds_df = data.frame(modelo1=preds1,modelo2=preds2,modelo3=preds3, real=real)
return(preds_df)
}## [1] "check this problem again"
## [1] "check this problem again"
## [1] "check this problem again"
## [1] "check this problem again"##
## Attaching package: 'plotly'## The following object is masked from 'package:MASS':
##
## select## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layout
Como podemos ver, los modelos 1 y 2 se “pisan” prácticamente, clavando el uno las predicciones del otro. Estos dos tienden a estimar algo por encima del valor real, como vemos en el gráfico de arriba, mientras que el tercer modelo suele estimar algo por debajo. Visualmente parece que el modelo que mejor predice es el 3; vamos a constatarlo en números.
Estadísticas de Error del Modelo 1
## ME RMSE MAE MPE MAPE
## Test set -809.9391 1527.874 1208.677 -11.54855 15.44051Estadísticas de Error del Modelo 2
## ME RMSE MAE MPE MAPE
## Test set -790.9721 1463.625 1153.604 -11.20356 14.70737Estadísticas de Error del Modelo 3
## ME RMSE MAE MPE MAPE
## Test set 297.8611 1222.87 996.9501 2.331691 11.63035Como decíamos, el modelo 3 es el que mejor predice de estos 3, ya que tiene un enor MAPE y menor RMSE (el MAPE nos permite medir el error en términos porcentuales, en este caso un 11.63% (está multiplicado por 100)); mientras que el RMSE nos permite saber cuánto nos equivocamos en las unidades que medimos, en este caso en el número de espectadores al mes en los cines. En el Mean Error vemos que los primeros dos modelos suelen quedarse por encima del valor real, pues este es negativo, mientras que el del tercero es positivo, mostrando que suele quedarse por debajo del valor real.
AUTOMATIC SARIMAX-GARCH: MACHINE LEARNING APPROACH TOWARDS TIME SERIES PREDICTIONS
En este pequeño apartado vamos a contrastar los resutados que obtuvimos anteriormente con el ajuste automático, al que además le vamos a incluir una modelización de la varianza, para estimar los valores superiores e inferiores posibles.
SARIMAX_GARCH = function(d=datos,nendtrain=96, var="users", datecol="date"){
preds = c()
mean_val = c()
results1 = c()
results2 = c()
for(i in seq(nendtrain, nrow(d)-12,12)) {
train = data.frame(d[1:i, ])
test = data.frame(d[(i+1):(i+12), ])
train.ts = ts(train[,var], start=c(year(train[1, datecol]), month(train[1, datecol])),
end=c(year(train[nrow(train), datecol]), month(train[nrow(train), datecol])),
frequency=12)
test.ts = ts(test[, var], start=c(year(test[1, datecol]), month(test[1,datecol])),
end = c(year(test[nrow(test), datecol]), month(test[nrow(test),datecol])),
frequency=12)
xregs = as.matrix(calculoExplicativasCalendario(train$date,0)[,c("semanaSanta", "dt", "bisiesto", "festivo", "price")])
ajusteAutomatico<-auto.arima(train.ts,max.d=1, max.D=1, max.p=1, max.q=1, max.P=1, max.Q=1,
seasonal=TRUE,ic="aic",xreg=xregs,stepwise=TRUE)
listaOutliers <- locate.outliers(ajusteAutomatico$residuals,
pars = coefs2poly(ajusteAutomatico),
types = c("AO", "LS", "TC"),cval=3)
outliers <- outliers(c(as.character(listaOutliers$type)), c(listaOutliers$ind))
outliersVariables <- outliers.effects(outliers, length(ajusteAutomatico$residuals))
RegressorsMasOutliers <- cbind(xregs,outliersVariables)
newxregs=as.matrix(calculoExplicativasCalendario(test$date, 0)[, c("semanaSanta", "dt", "bisiesto", "festivo", "price")])
outliersNew <- outliersVariables[(i-11):i,]
newxregs <- as.matrix(cbind(newxregs,outliersNew))
#test = btc2$btc[i:i+5] #we set this to have 5 points so that the out.sample can work.
m1.mean.model <- auto.arima(train.ts,max.d=1, max.D=1, max.p=1, max.q=1, max.P=1, max.Q=1,
seasonal=TRUE,ic="aic",xreg=RegressorsMasOutliers,stepwise=TRUE)
predicc = predict(m1.mean.model, n.ahead=12, newxreg= newxregs)
preds = c(preds, predicc$pred)
ar.comp <- arimaorder(m1.mean.model)[1]
d_ <- arimaorder(m1.mean.model)[2]
ma.comp <- arimaorder(m1.mean.model)[3]
library(rugarch)
model.garch = ugarchspec(mean.model=list(armaOrder=c(ar.comp,d_, ma.comp), arfima =F),
variance.model=list(garchOrder=c(1,1)),
distribution.model = "std")
model.garch.fit = ugarchfit(data = train.ts, spec=model.garch, solver = 'lbfgs' )
modelfor=ugarchforecast(model.garch.fit, data = NULL, n.ahead = 12)
mean_val = c(mean_val, modelfor@forecast$seriesFor[,1])
results1 <- c(results1, predicc$pred + modelfor@forecast$sigmaFor[,1])
results2 <- c(results2, predicc$pred - modelfor@forecast$sigmaFor[,1])
}
real = data.frame(d)[(nendtrain+1):nrow(d), var]
df = data.frame(pred=preds,meanval=mean_val,upper=results1, lower=results2, real=real)
return(df)
}
Medidas de Error para el Modelo Automático
## ME RMSE MAE MPE MAPE
## Test set -729.9265 1650.171 1321.42 -11.12311 17.00692Vemos que el ajuste automático deja unos resultados bastante pobres en comparación con el mejor modelo de los anteriores, que es el que hacemos con \(y^\lambda\). En aquél caso teníamos un 11.65 de MAPE, aquí tenemos un 17… Por lo tanto nos quedaríamos con el modelo 3 del apartado anterior, pues éste es el que mejor predice. Vamos a volver a dibujar el gráfico de las predicciones de una forma más completa, introduciendo también las predicciones del modelo automático.